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Abstract We study the behavior of orbits in two different galactic dynamical models, 
describing the motion in the central parts of a triaxial elliptical galaxy with a dense nu- 
cleus. Numerical experiments show that both models display regular motion together with 
extended chaotic regions. A detailed investigation of the properties of motion is made for 
the 2D and 3D Hamiltonian systems, using a number of different dynamical parameters, 
such as the Poincare surface of a section, the maximal Lyapunov Characteristic Exponent, 
the S(c) spectrum, the S(w) spectrum and the P(f) indicator. The numerical calculations 
suggest that the properties of motion in both potentials are very similar. Our results show 
that one may use different kinds of gravitational potentials in order to describe the motion 
in triaxial galaxies while obtaining quantitatively similar results. 
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1 INTRODUCTION 



About forty years ago, the prevailing view was that elliptical galaxies were oblate spheroids flattened 
by rotation (see Sandage et al.|1970 1. It was the pioneering work of Bertola & Capaccioli| ( |1975[ ) and 
Illingworth ( 1977| l that led astronomers to abandon the assumption that elliptical galaxies are necessarily 
oblate. 

It is well known that observations of elliptical galaxies yield only the projected isophotes and, thus, 



determination of their intrinsic shapes requires statistical analysis, based on large samples (see Ryden 
1996 Alam & Ryden|2002 Vincent & Ryden|2005|l or mapping of potentials via detailed kinematical 
data for individual galaxies (see |Davies et al.|2001| Rest et al.|200T||Statler et al.|2004) . 

Today it is believed that the shapes of elliptical galaxies are prolate or triaxial rather than oblate 
(see |Benacchio & Galletta|1980||Binggeli|1980{|Alam & Ryden|2002] l. On the other hand, kinematical 
studies of elliptical galaxies show evidence that triaxial galaxies do exist. Moreover, observational data 
indicate that most of the triaxial elliptical galaxies host a black hole or a dense nucleus in their centers 
(see |Bak & Sta tler 2000. [Staffer et al!|2004) . On this basis, we believe that it would be of interest to 
investigate the dynamical properties of a triaxial elliptical galaxy, particularly in its central region. 

In order to describe the motion in the triaxial elliptical galaxy we use the well known logarithmic 
potential 
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where vq is used for the consistency of the galactic units, a and b are flattening parameters, and Cb is the 
scale length of the bulge component (see Binney & Tremaine 2008 ). Expanding potential (1) in a Taylor 
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series about the origin and keeping terms up to the fourth degree in the variables, we find 



K = 



lnc? + 



2c 2 
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Thus the polynomial potential V\ is 



where it was assumed that 
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The reader can find more details about the Taylor expansion of the logarithmic potential in Zotos 
( |201 lb| l. To potentials (1) and (3) we add the potential of a spherically symmetric nucleus 



V„ = 



-Af„ 



^/x 2 + y 2 + z 2 + c 2 ' 



(5) 



where M n is the nuclear mass, while c n is the scale length of the nucleus. We apply a Plummer sphere, in 
order to increase the central mass of the galaxy. This method has been applied several times in previous 
works, having as an objective to study the effects of the introduction of a central mass component in a 
galaxy (see Hasan & Norman|1 990; Has an et al.|1993[ >. 

The aim of the present article is to investigate the properties of motion near the center of a triaxial 
elliptical galaxy described by potentials (1) and (3) with an additional dense nucleus described by po- 
tential (5). In particular, we are interested to study the regular or chaotic character of motion in both of 
the above described potentials, in order to be able to compare the corresponding results. Furthermore, 
we shall compare the density in the central parts of the triaxial galaxy derived using the two potentials 
described above. In order to achieve a better picture for the properties of motion, we first investigate the 
2D system, that is when z — 0, and then we will use the corresponding results to study the dynamical 
system of three (3D) degrees of freedom. 



From the pioneering work of |Henon & Heiles ( 1964 1 there has been an ongoing interest in finding 
new methods, in order to distinguish between ordered and chaotic motion in dynamical systems. The 
Poincare surface of section (PSS) for the two dimensional (2D) systems and the Lyapunov Characteristic 
Exponent (LCE) ( Benettin et al.|1976 Froeschle|1984| Lichtenberg & Lieberman|1992| l for dynamical 
systems with any degree of freedom are two well known methods to characterize an orbit as regular or 
chaotic. Over the last thirty years, an effort has been made in order to find new, modern and also reliable 
and fast ways to detect the chaotic behavior in galactic systems. One could mention the frequency map 
analysis developed by Laskar ( Laskar et al.||1992) Laskar 1993), the dynamical spectra of stretching 
numbers (the distribution of values of a given parameter along the obit), introduced by Froeschle ( 1984) 
(see also |Froeschle et al.|[T993| |Voglis & Contopoulos||1994 |Contopoulos et aT1|1995[ |Contopoulos| 
et al.|1997| l and the P(f) spectral method applied and used by Karanis & Vozikis (2008). In the present 
research, we use, apart from the classical PSS technique and the LCE, some modern methods such as 
the S(c) and S(w) dynamical spectra and the P(f) indicator. 

Here, we must provide some additional theoretical information regarding these new dynamical 
methods. We use the S(c) spectrum in order to characterize the nature of an orbit in a 2D dynami- 



cal system. This spectrum has been proved to be a very reliable tool in several cases (see Caranicolas & 



Papad opoulos|2007 Zotos|20PTa I. The nature of a 2D orbit can be revealed by looking at the shape of 
the S(c) spectrum. If the shape of the spectrum is a well defined [/-type structure, then the correspond- 
ing orbit is regular. On the other hand, if the shape is complicated and highly asymmetric, with a lot of 
large and small abrupt peaks, then the orbit is chaotic. Moreover, the S(c) spectrum can help us identify 
resonant orbits of higher multiplicity, as it produces as many [/-type structures as the total number of 
islands of the invariant curves on the x — p x phase plane. One more advantage of this spectrum is that it 
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can be deployed in order to calculate the sticky period of a 2D orbit and also to follow its time evolution 
towards the chaotic sea (see fig. 5 in |Zotos|2011a| ). 

For the study of 3D orbits, we use the S(w) spectrum. By definition this spectrum is based on a 
complicated combination of the coordinates and the momenta of the 3D orbit. In particular, this spectrum 
is an advanced form of the 5(c) spectrum and therefore carries all the characteristics mentioned in 
the previous paragraph regarding the pattern of the spectrum for regular and chaotic orbits. The only 
difference is that in this case the S(w) spectrum produces as many U -type structures as the total number 



of invariant 3D tori in the (x,p x , z) phase space. We introduced this new spectrum definition in Zotos 
(201 la|l, in order to construct a new spectral definition appropriate for the study of 3D orbits 



The Fourier Transform is usually defined as a transformation of a quantity q which is a function 
of time, q(t), with regard to its respective function of amplitude p, which is a function of frequency 
P(f). We can define a series of time intervals between successive crossings over a section. Then we 
calculate the Power Spectrum of these time intervals, using a Discrete Fast Fourier Transformation 
(FFT) algorithm. Looking at the P(f) spectrum of a regular 2D or 3D orbit, we expect to observe a 
smooth curve, with some additional peaks corresponding to the "periodicities" of the time series. On the 
contrary, in a 2D or 3D chaotic orbit no such "periodicities" exist and therefore its P(f) spectrum will 
produce a very "noisy" pattern with a large number of peaks which would be very densely distributed. 
One of the main advantages of this spectral method is that is uses only one orbit and we do not need 
to trace the behavior of any nearby orbit. Furthermore, the detection can be made quite early using less 
iterations, compared to the iterations needed to reach a conclusive result using the LCE. More detailed 
information regarding this method and its applications can be found in Karanis & Vozikis (2008 ). 



Here we must remind the reader that the 5(c) spectrum is the distribution function of the parameter 

where AN(c) are the numbers of the parameters c in the interval (c, c + Ac) after N iterations. The 
parameter c is defined as 

X{ Pxi ,~\ 

Ci = , (7) 

Pyi 

where (xi,p X i,p y i) are the successive values of the (x,p x ,p y ) elements of the 2D orbits, on the Poincare 
x — p x , y = 0,p y > phase plane. More details regarding the 5(c) spectrum and its applications can 



be found in Caranicolas & Papadopoulos (2007) and Caranicolas & Zotos (2010). 



The Hamiltonian corresponding to the potential (1) or (3) is written as 

H=^(p 2 x +pl+p 2 z )+V t (x,y,z) = E, (8) 

where V t represents V tg = V g + V n or V t i = V\ + V n . Here p x ,p y and p z are the momenta per unit 
mass conjugate to x, y and z respectively, while E is the numerical value of the Hamiltonian (8), which 
is conserved. 

In this article, we use a system of galactic units, where the unit of length is lkpc, the unit of mass 
is 2.325 x 10 7 Af© and the unit of time is 0.97748 x 10 s yr. The velocity unit is 10 km s -1 , while G is 
equal to unity. In the above units we use the values: vo — 10, Cb = 3, M n — 10, c n = 0.1, a = 1.5 and 
6 = 1.7. 

The results of the present research are based on the numerical integration of the equations of motion 

.. _ dV t (x,y,z) 
dx 

dV t (x,y,z) 

y = a ' 

ay 

z = _ g%M , (9) 
oz 
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which was made using a Bulirsh-Stoer routine in Fortran 95, with double precision in all subroutines. 
The accuracy of the calculations was checked by the constancy of the energy integral (8), which was 
conserved up to the eighteenth significant figure. 

This paper is organized as follows. In Section 2 we present and compare the results for the 2D 
systems. In Section 3 we compare the mass density near the center derived using the two 3D potentials. 
Moreover, we compare the properties of the 3D orbits in the two dynamical models. In Section 4, the 
conclusions and the discussion of our results are presented. 

2 RESULTS FOR THE 2D DYNAMICAL SYSTEMS 

In this section, we study the character of orbits in the 2D dynamical systems. In this case, we set z = 
p z = in (8) and the corresponding 2D Hamiltonian is 

H 2 = \ (pI + P 2 y ) + V t (x, y) = E 2 , (10) 

where E% is the numerical value of the Hamiltonian. As the phase space of the system is four dimen- 
sional we use the x — p x ,y — 0,p y > Poincare phase plane. The results are presented in Figures [T] 
and|2 

Figure [T] shows the phase plane for potential V tg , while Figure [2] shows the phase plane for poten- 
tial Vti. Here, we must emphasize that the two phase planes were constructed for values of energies 
connected by the relation 

E 2tg = E 2t i + vl\nc h (11) 

and the same initial conditions, in order to be able to make the comparison. Here we took E 2t \ — 
—4.70, Vq lncb = 109.86, which gives the value 105.16 for E 2tg . As one can see, the two phase planes 
are almost identical. In both figures, we see areas of regular motion and extended chaotic regions. There 
are three main families of regular orbits, (i) Orbits producing invariant curves around each of the two 
stable 1:1 resonant periodic points, (ii) Orbits producing a set of three islands - one of them is on the 
x-axis, while the other two are symmetric with respect to the p x axis. These orbits are characteristic of 
the 3:3 resonance, (iii) Box orbits producing invariant curves surrounding the whole chaotic sea. Note 
that the area on the x — p x phase planes occupied by each of the above families of regular orbits is 
quantitatively the same in both Figures [T] and [2] In addition to the regular orbits there are also a large 
number of irregular orbits producing a large, unified chaotic sea. Note that the extent of the chaotic sea 
is almost the same in both phase planes. The differences between the two phase planes produced by 
potentials Vt g and Vti are negligible and they are confined to some tiny islands, embedded inside the 
chaotic sea. These tiny islands are produced by secondary resonances. 

In order to investigate and compare in detail the properties of motion in both potentials, we present 
and compare in the following a number of orbits belonging to different families of orbits. 

Figure [3|a)-(d) shows results for a regular orbit in potential Vtg. The orbit shown in Figure |3j a) 
belongs to family (i) and has initial conditions: xq — 0.5, yo — 0,p x o = 0, while in all cases, p y o is 
found from the energy integral (10). The corresponding values of energy and all the other parameters 
are as in Figure [T| Figure[3|b) shows the maximum LCE of the orbit, which vanishes indicating regular 
motion. FigurepFc) shows the S(c) spectrum of the orbit. Here, we see a well defined U type spectrum 
characteristic of the regular motion. In Figure |3jd), we see a plot of the P(f) indicator, which displays 
only two peaks indicating regular motion. In order to help the reader, we note that the orbit shown in 
Figure[3ja) was calculated for a time period of 100 time units. The time scale for the S(c) spectrum and 
the P(f) indicator was 10 3 — 10 4 time units. 

Figure |4ja)— (d) shows results for an orbit, with the same initial conditions and with the same time 
scales for all calculations, but for the potential Vti. All other parameters are the same as in Figure [2] 

In order to have a better picture of the properties of motion in the two potentials Vt g and Vti, we 
present results for two more orbits. 

Figures |5ja)-(d) and |6ja)— (d) are similar to |3ja)-(d) and |4f a)— (d), for an orbit with initial condi- 
tions: x n = 0.15, j/o = ®,PxO = 4.5. This orbit belongs to family (ii) and it is characteristic of the 
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3:3 resonance. As one can see, the outcomes presented in the two figures are very similar. The results 
presented in Figures |7ja)— (d) and[HJa)-(d) are similar to those of Figures |3ja)-(d) and|4|a)-(d) but for a 
chaotic orbit. Initial conditions are: xo = 0.02, yo = 0,p x a = 2.5. Here the maximum LCE has a pos- 
itive value indicating chaotic motion. Moreover, the S(c) spectrum shows a number of large and small 
peaks which is characteristic of chaotic motion. Finally, the P(f) indicator is highly asymmetric with a 
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large number of peaks, also indicating chaotic motion. Comparing the results given in Figures |7Ja)-(d) 
and[SJa)-(d), we can say that they are very similar. 

Given all the above, we can say that our numerical results, which are obtained by several different 
dynamical methods using regular and chaotic orbits, strongly suggest that the potential Vti satisfactorily 
describes the properties of motion of the potential Vt g near the center of a triaxial, elliptical galaxy. 
Since the potential of the spherical nucleus V n is the same in both potentials V tg and Vti, this means that 
no information is lost when we go from a global triaxial logarithmic potential (1) to the local polynomial 
potential (2). Remember that this only holds near the center of the galaxy, when (4) is valid. In order to 
investigate and compare the character of motion in the two 2D potentials Vt g and Vti, we have computed 
a large number of orbits - about 1000 - with different initial conditions (xq, p x o), but with the same initial 
conditions in both 2D potentials. In particular, as we have in both cases regular regions and only one 
unified chaotic sea in each x — p x phase plane, we calculate the maximum value of the LCE by choosing 
500 orbits with different and random initial conditions (xo,p x o) in the regular regions and 500 orbits 
with different and random initial conditions (x ,p x0 ) in the chaotic sea in each case. Our numerical 
experiments show that the vast majority of orbits - about 97.4% - displayed the same characteristics, 
including the same nature of orbit, maximum LCE, S(c) spectrum and P(f) indicator, while only 2.6% 
of the tested orbits were different. 



3 RESULTS FOR THE 3D DYNAMICAL SYSTEMS 

We now proceed to study the properties of motion in the 3D potentials. Before doing this, it would be 
interesting to compare the mass densities derived from the 3D potentials Vtg and Vti. The mass density 
can be found using Poisson's law 

V 2 14 = AnGpt , (12) 

where Vt represents Vt g — V g + V n or Vti = Vi + V u , while p t represents p tg or p t \. 

Figure |9|a)-(h) shows the surfaces of equal density for the 3D potentials Vt g and Vti- We can see 
that the results are very similar. In order to compare the mass density from another point of view, we 



present in Figure 10 (a)-(f) the contours of equal density in the xy, xz and yz planes respectively for the 
two potentials. As is understood, these contours are the projections of the four surfaces of equal density 
to the three principal planes xy, xz and yz. Here we can visualize that the deviations between the mass 
density of the two potentials are extremely small and therefore negligible. 

Let us now come to investigate and compare the orbits in the two 3D potentials. For this purpose, 



we apply the S(w) dynamical spectrum, which was introduced in Zotos (201 la i, in order to distinguish 



between ordered and chaotic motion in 3D dynamical systems. The parameter W{ is defined as 

(Xi - p xi ) - (Zi - p zi ) 

Wi = , (13) 

Pyi 

where [xi, Zi,p X i,p y i,p Z i) are the successive values of the (x, z,p x ,p y ,p z ) elements of the 3D orbits. 
The dynamical spectrum of the parameter w is its distribution function 

NAw 

where AN(w) is the number of parameters w in the interval (w, w + Aw) after N iterations. In order 
to study the character of a 3D orbit, the S(c) spectrum can also be applied. Note that the coupling of the 
third component z, carrying all the information regarding the 3D motion, is hidden in the definition of 
the S(c) spectrum, but in any case it affects the values of x, p x and p y . Using the definition of the S(w) 
spectrum, we overcome this minor drawback as we deploy an improved dynamical spectrum, especially 
suitable for 3D orbits. 

Figur e [TT| a)-(d) shows the results for a 3D regular orbit in potential Vt g . The orbit which is shown 
in Figure [Tlfa), has initial conditions: xq — 0.5, yo = p x o = p z o = 0, z$ = 0.1, while for all 3D 
orbits the value of p y0 is found from the energy integral (8). The corresponding values of all the other 
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parameters are the same as in Figure [T] The value of energy is E tg = 105.16, the same as in the 2D 
system. The maximum LCE of this orbit, which is shown in Figure [TTJb), vanishes indicating regular 
motion. Figure 1 1 'c) shows the S(w) spectrum of the orbit. This is a well defined U type spectrum 
characteristic of the regular motion. In Figure [TT^d) we can see the P(f) indicator which also indicates 
regular motion. Figure [T2|a)-(d) shows results for the same orbit but in the potential Vti- The values of 
the other parameters are the same as in Figure[2] The value of energy is E t i = —4.70, the same as in the 
2D system. Comparing the two Figures [TT|a)-(d) and[T2|a)-(d) we see that the results are very similar. 

Figures 13 a)-(d) and 14 'a)-(d) are similar to Figures [TTJa)-(d) and 12 a)-(d) for a resonant 3D 
orbit with initial conditions: xq = 0.15, yo = 0,p x o = 4.5, p Z Q — 0, zq — 0.01. The similarity between 



the two patterns is evident. Finally, in Figures 15 a)-(d) and[T6|a)-(d), we present results for a chaotic 
3D orbit. The initial conditions are: xq = 0.02, yo = 0,^0 = 2.5, p z o — 0,zq — 0.1. The values of 



energy and other parameters are as in Figure 13 a)-(d) and 14 a)-(d) respectively. Again we see that the 
results are very similar. 

In order to investigate and compare the characteristics of motion in the two 3D potentials V tg and 
Vti, we work as follows. We use initial conditions (xo, p x o, zq), yo = p z o = 0, where (xrj, p x o) is a point 
on the phase planes of the corresponding 2D potentials. This point lies inside the limiting curve, which 
is the curve containing all the invariant curves of the 2D system. The equation of the limiting curve is 



1 



p 2 x + V t {x) = E 2 



(15) 



where E% is i?2tg or -E^tl- Using this method, we have computed a large number of 3D orbits - about 1000 
- with the same initial conditions in both 3D potentials Vtg and Vti. In particular, as we have in both cases 
regular regions and only one unified chaotic sea in each x — p x phase plane, we calculate the maximum 
value of the LCE by choosing 500 orbits with different and random initial conditions (xo,p x o, zo) in the 
regular regions and 500 orbits with different and random initial conditions (xo,p x o, z ) in the chaotic 
sea in each case. Our numerical calculations indicate that the majority of orbits - about 94.6% - displayed 
almost the same characteristics, which are the shape of the orbit, the maximum LCE, the S(w) spectrum 
and the P(f) indicator, while only 5.4% of orbits were different. 

Therefore, from the investigation of the 3D potentials, we have arrived at the following conclusions. 
The mass densities near the center of the elliptical galaxy produced by the potentials are nearly the same. 
Furthermore, orbits with the same initial conditions in both potentials are very similar and show similar 
patterns of the maximum LCE, the S(w) spectrum and the P(f) indicator. Moreover, the percentage 
of chaotic orbits in both 3D potentials seems to be almost the same. Thus, we conclude that, generally 
speaking, the properties of motion in both 3D potentials are almost the same. 



4 DISCUSSION AND CONCLUSIONS 

In this paper we have studied the properties of motion near the center of a triaxial elliptical galaxy 
described by two different potentials Vt g = V g + V n or Vti = Vi + V n . In fact, V\ is an expansion of the 
potential V g in a Taylor series near the center, up to the terms of fourth degree in the variables, while 
the potential V n was added for two basic reasons. The first reason is that there is observational evidence 
that black holes or dense massive nuclei lie in the centers of some elliptical galaxies. The second reason 
is that with the additional term V n , potentials Vt g and Vti produce interesting orbital characteristics, 
such as several families of periodic orbits together with large chaotic regions. In this work we do not 
have as an objective to provide anything new regarding the properties of motion of these dynamical 
systems. On the contrary, we use well known potentials and try to compare them by using different 
kinds of indicators. Our purpose is to show how we can correctly expand a logarithmic potential in a 
Taylor series to produce a harmonic oscillator. The main result from our research is that despite the fact 
that the potentials are different, they display almost identical properties of motion. The results obtained 
using different dynamical indicators are very similar. This means that the Taylor expansion is valid 
and the harmonic oscillator potential can satisfactorily describe the local motion in the central parts of 
an elliptical galaxy. On this basis, we provide relations regarding the involved parameters, so that the 
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parameters of the system do not have arbitrary values but rather values which are related to the global 
logarithmic potential and they have physical meaning. 

First we studied the 2D system. The phase planes which were constructed for the two above dif- 
ferent potentials were found to be nearly identical. In the next step we studied the properties of orbits 
with the same initial conditions in both potentials using the maximum LCE, the 5(c) spectrum and the 
P(f) indicator. In all cases the results were very similar. Then we started the study of the 3D system 
by comparing the mass density in the two potentials Vt g and V%\. The results have shown very small 
differences in the mass densities. As in the 2D system we also investigated the properties of orbits in 
both 3D potentials using the maximum LCE, the S(w) spectrum and the P(f) indicator. The results 
were once more very similar. Furthermore, our numerical calculations suggest that the percentage of 
chaotic orbits is about the same in both the above potentials. 

Also note that, strictly speaking, potential (1) is a global galactic potential, which describes a triaxial 
galaxy as a whole, while potential (2) is a local potential, which describes the galaxy only in its central 
parts. In other words, the description is satisfactory only if relation (4) is valid. Since the two potentials 
are similar, the orbital behavior of the orbits should be almost identical, while the minor observed 
differences are caused by the higher order terms of the Taylor expansion. It would be of particular 
interest to inspect and locate the range of the parameters for which the orbital behavior in both dynamical 
systems remains the same. Numerical experiments indicate that the results are sensitive to the parameters 
of the dynamical systems. In particular, we conclude that the properties of motion (2D or 3D) in both 
potentials V tg and Vti are almost the same only when 1.1 < a < 1.9,0.1 < b < 1.8,1.8 < Cb < 
3.2,5 < M D < 25 and 0.10 < c n < 0.25. Numerical calculations not given here show that the 
properties of motion near the center of potentials (1) and (3) are almost the same. The only difference is 
that in this case (when the spherical nucleus in not present) we only observe regular motion, while the 
chaotic orbits if any are negligible. With the additional term V n , the two potentials display regular and 
chaotic motion as well and the properties of motion are again very similar. 

Here we must remind the reader that he or she can find the definitions and also some useful the- 



oretical explanations about the 5(c) and S(w) dynamical spectrums in Caranicolas & Papadopoulos 



( |2007| l, in |Caranicolas & Zoto s (2010) and also in Z otos| (|201 la[>. The defi nition and additional infor- 
mation regarding the P(f) indicator are given in |Karanis & Vozikis| ( |2008) l. The main drawback of all 
these methods is that they can only provide qualitative results regarding the regular or chaotic nature of 
an orbit. Therefore, we must check the shape of the indicator by eye each time in order to characterize 
an orbit. Nevertheless, these dynamical indicators are very useful as they can provide fast and reliable 
results. In order to check their validity and reliability in each case (2D and 3D systems), we have com- 
pared these qualitative results with a highly accurate and quantitative method, such as the Lyapunov 
Characteristic Exponent. Our comparison proves that although the outcomes of these spectral methods 
are qualitative they are also very reliable. 
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Fig. 3 (a) A regular orbit in the 2D potential V tg . Initial conditions are: x = 0.5, yo = 
0,p X Q = 0, while p y o is found from the energy integral. The values of all other parameters 
and energy are as in Fig.[T] (b) A plot of the maximum LCE vs. time for the orbit shown in 
(a), (c) The S(c) spectrum of the orbit shown in (a) and (d) The P(f) indicator for the orbit 
shown in (a). 
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Fig. 4 (a)-(d): Similar to Fig.pja)-(d) for the potential V%\. The values of all other parameters 
and energy are the same as in rig. [2] 
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Fig. 6 (a)-(d): Similar to Fig.pfa)-(d) for the potential Vy. The values of all other parameters 
and energy are the same as in rig. [2] 
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Fig. 9 (a)-(h): Surfaces of equal density for the 3D potentials. Left patterns correspond to 
potential V tg , while right patterns to potential V t \. 
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Fig. 10 (a)-(f): Contours of equal density in the xy, xz and yz planes. Left patterns 
spond to potential V tg , while right patterns to potential V t \. 
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Fig. 11 (a) A regular orbit in the 3D potential Vt g . Initial conditions are: xq = 0.5, ya = 
®,Px0 = 0, z Q = 0.1, while p y0 is found from the energy integral. The values of all other 
parameters and energy are as in Fig. [T] (b) A plot of the maximum LCE vs time for the orbit 
shown in (a), (c) The S(w) spectrum of the orbit shown in (a) and (d) The P(f) indicator for 
the orbit shown in (a). 
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Fig. 13 (a)-(d): Similar to Fig. 11 'a)-(d) for a resonant 3D orbit. Initial conditions are: xq = 
0.15,yo = 0,p x o = 4.5, zo = 0.01. The values of all other parameters and energy are the 
same as in Fig.[T] 
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Fig. 15 (a)-(d): Similar to Fig. 1 1 'a)-(d) for a chaotic 3D orbit. Initial conditions are: x a = 
0.02, yo — 0, p x o = 2.5, zo = OTTTThe values of all other parameters and energy are the same 
as in Fig.[T] 



402 



Euaggelos E. Zotos 




